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Abstract 

We investigate the quantum chaotic properties of the Dicke Hamiltonian; a quantum-optical 
model which describes a single-mode bosonic field interacting with an ensemble of N two-level 
atoms. This model exhibits a zero-temperature quantum phase transition in the N — > oo limit, 
which we describe exactly in an effective Hamiltonian approach. We then numerically investigate 
the system at finite N and, by analysing the level statistics, we demonstrate that the system under- 
goes a transition from quasi-integrability to quantum chaotic, and that this transition is caused by 
the precursors of the quantum phase-transition. Our considerations of the wavefunction indicate 
that this is connected with a delocalisation of the system and the emergence of macroscopic coher- 
ence. We also derive a semi-classical Dicke model, which exhibits analogues of all the important 
features of the quantum model, such as the phase transition and the concurrent onset of chaos. 
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I. INTRODUCTION 



Chaos plays a key role in considerations concerning the boundary between the classical 
and quantum worlds, not just because of the importance of chaos in classical physics jlj], 
but becaus there is no direct analogue of chaos in quantum mechanics The linearity of 
quantum dynamics precludes the characteristic exponential sensitivity to initial conditions of 
classical chaos, and forces us to look for what have become known as "signatures of quantum 
chaos" - properties whose presence in the quantum system would lead us to expect the 
corresponding classical motion to be chaotic jal . Several such signatures have been identified, 
such as level statistics jj, 0], level dynamics p, and sensitivity to initial perturbation 

An oft encountered feature of quantum-chaotic systems is that as some parameter is 
varied, these signatures bespeak a cross-over from integrable to quantum chaotic behaviour. 
This parameter may, for example, describe the character of boundary conditions, such as the 
shap e of a q uantum billiard [5(], the distribution of random fluctuations in disorder models 

i3Q,Q,Q,Qli2i. 



, or the strength of some non-linear potential or interaction 
A large class of model may be described by a Hamiltonian of the form 



H = H + XV, (1) 



where, although Hq is integrable, the full Hamiltonian H is not for any A 7^ 0. Here, increas- 
ing the parameter A from zero upwards gradually drives the system away from integrability 
and towards chaos. A well studied, albeit time-dependent, example is the kicked rotator |3[, 
where the parameter A is the kick-strength. 

In this paper, we consider a system of the type described by Hamiltonian but un- 
like the typically one-dimensional or non-interacting models, we shall consider a system of 
iV interacting particles, in a situation where many-body and collective effects are critical. 
Specifically, the model we study exhibits a quantum phase transition (i.e. one at zero tem- 
perature 18[) in the thermodynamic limit of N — ► 00 at a critical value of the parameter, 
A c . 

The influence of a quantum phase transition (QPT) on the transition to chaos has been 
studied in but a handful of cases. Important examples include the three-dimensional An- 
derson model, where the metal-insulator transition is accompanied by a change in the level- 
statistics jsj], and models of spin glass shards |lfl| . which have found topical application in 
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the study of the effects of quantum chaos on quantum computing If] . Heiss and co-workers 
have investigated the connection between the onset of chaos near a QPT and the exceptional 
joints of the spectrum , both generically and for the specific example of the Lipkin model 



In order to investigate the impact of QPT on the signatures of quantum chaos, we study 



the Dicke Hamiltonian (DH) 21 , which is of key importance as a model describing collective 
effects in quantum optics . We demonstrate that there is a clear connection between 

the precursors of the QPT and the onset of quantum chaos as manifested in the level- 
statistics. We are able to understand this connection by studying the wavefunctions of the 
system, and by deriving a semi-classical analogue of this intrinsically quantum system. The 
current publication is an extension of our previous work [24J]. 

In the form considered here, the DH describes a collection of N two-level atoms interacting 
with a single bosonic mode via a dipole interaction with an atom-field coupling strength A. 
The DH may be written 

H = Huj J z + huja^a + — L= (a ] + a) (J. + J_) , (2) 

V 2 J 

where a, describe a bosonic mode of frequency ou, and the angular momentum operators 
{ Jf, i = z, ±} describe the ensemble of two-level atoms of level-splitting ujq in terms of a 
pseudo-spin of length j = N/2. The thermodynamic limit of iV — > 00 is thus equivalent to 
making the length of the pseudo-spin tend to infinity j — > 00. The DH is usually considered 
in the standard quantum optics approach of the rotating- wave approximation (RWA), which 
is valid for small values of the coupling A, and involves neglecting the counter-rotating terms 
a) J + and aJ_. This makes the DH integrable, simplifying the analysis but also removing 
the possibility of quantum chaos. Dicke used this model to illustrate the importance of 
collective effects in the atom-light interaction |2^], leading to the concept of super-radiance, 
where the atomic ensemble spontaneously emits with an intensity proportional to N 2 rather 
than N, as one would expect if the atoms were radiating incoherently j^ . 

The phase transition in the DH was first described by Hepp and Lieb J25| , and a math- 
ematically more transparent treatment was provided by Wang and Hioe |26|. They consid- 
ered the thermodynamics of the model in the RWA and concluded that for a coupling of 
A < ^ujuq, no phase transition occurs for any temperature, whereas for A > ^ujujq, there 
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exists a critical temperature T c given by 

kk = ^ artanh (99 ' (3) 

at which point the system undergoes a phase transition. Above the critical temperature, 
the system is in the effectively unexcited "normal phase", whereas for T < T c the system 
is in the "super-radiant phase" , a macroscopically excited and highly collective state which 
possesses the potential to super-radiate. 

In contrast to this earlier work, we shall consider this phase transition at zero temperature, 
where increasing the coupling A through a critical value of A c = drives the system to 

undergo a transition from the normal to the super-radiant phase (the difference between this 
critical coupling A c and the value quoted for the finite-temperature case arises because the 
latter has been derived in the RWA, which renormalises the critical coupling by a factor of 



two 
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281]). Here, we derive exact results without the RWA for the energy spectrum and 



eigenfunctions in the thermodynamic limit by employing a bosonisation technique 



Dased 
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upon the Holstein-Primakoff transformation of the angular momentum algebra 
This enables us to derive an effective Hamiltonian to describe the system in each of its 
two phases. One important step that we make is the introduction of an abstract position- 
momentum representation for both the field and atomic systems. This not only facilitates the 
formulation of the exact solutions, but also provides us with a useful way of visualising the 
wavefunctions across the phase transition. There is a discrete "parity" symmetry associated 
with this model, and at the phase-transition this symmetry becomes broken. This QPT has 



been discussed in the RWA by Hillery and Mlodinow [3lj . using an effective Hamiltonian 
method that is similar to ours. However, having illustrated the existence of the QPT, they 
concentrated solely on the normal phase, and were not interested in chaos. 

Away from the thermodynamic limit at finite iV and j, the DH is, in general, non- 
integra ble. Quantum-chaotic properties of the DH have been discussed by several authors 
I34I I35I 3(1 37, 38, 3^ but, to the best of our knowledge, have never been connected 



3J,|3i 



with the QPT, and a systematic study of the dependence of the systems behaviour on the 
number of atoms N is lacking. Graham and Hohnerbach have contributed extensively to the 
discussion j^, especially in relation to the special case of spin-1/2 (the Rabi Hamiltonian), 
and have outlined many semi-classical and approximate schemes for these systems. More- 
over, they have provided a preliminary analysis of the level statistics of the DH, concluding 
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that spectra of the type associated with quantum chaos do occur for certain, isolated param- 
eter values j^J. Several authors have conducted studies of chaos in various (semi-)classical 



models related to DH 



36 



37 



|38|, 39]. different semi-classical 

models is a consequence of the ambiguity in describing quantum spins in classical terms. 
The influence of the QPT also seems to have have been overlooked in these semi-classical 
models. 

We consider the quantum-mechanical system away from the thermodynamic limit by 
using numerical diagonalisation, and examine the energy spectra of the system for signatures 
of quantum chaos. We consider the nearest-neighbour level-spacing distribution function 
P(S), which is perhaps the best-known signature of quantum chaos We calculate the 
P (S) for various values of iV and A and demonstrate a clear connection between the change 
in P (S) from quasi-integrable to quantum chaotic and the coupling at which the QPT 
occurs, A c . We then proceed to consider the wavefunctions of the system at finite N using an 
abstract position-momentum representation. This enables us to conclude that the precursors 
of the QPT give rise to a localisation-delocalisation transition in which the ground-state 
wavefunction bifurcates into a macroscopic superposition for any N < oo. 

As mentioned above, there has been much work in trying to find a semi-classical ana- 



logue of the DH 



36 



37 



The bosonisation procedure that we employ here allows 



us to write the DH in terms of a pair of coupled harmonic oscillators. This suggests a 
very natural semi-classical analogue of the DH, obtained by simply replacing the quantum 
oscillators with classical ones. We demonstrate that our semi-classical model reflects the 
quantum behaviour better than those of previous studies. Specifically, our semi-classical 
model exhibits a symmetry-breaking phase transition in the limit that iV — > oo, and we 
show that the precursors of this classical transition give rise to the onset of classical chaos, 
in close agreement with the quantum model. An analogue of the macroscopic superposition 
is also evident. In our conclusions, we pay special attention to the meaning of a classical 
limit for the DH, and in particular the relevance of the semi-classical model derived here. 

The paper is organised as follows. In section |H] we introduce the DH fully. Exact solutions 
are derived in the thermodynamic limit in section HTTl Section HVl sees an analysis of the level- 
statistics and wavefunctions of the system at finite j. Our semi-classical model is derived 
in section and its phase transition and chaotic properties discussed. We discuss briefly 
the differences between the full DH and the Hamiltonian in RWA in section IVI| before we 
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draw our final conclusions in section IVHI Some of our exact expressions are reproduced in 
the Appendix. 
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II. THE DICKE HAMILTONIAN 



The full Dicke Hamiltonian (DH) models the interaction of N atoms with a number of 
bosonic field modes via dipole interactions within an ideal cavity [21]. We initially represent 
the atoms as a collection of TV identical, but distinguishable two-level systems each with level- 
splitting co>o- The zth atom is described by the spin-half operators |sj^; k = z, ± j, obeying 
the commutation rules [s z , s±] = ±s±; [s+, s_] = 2s z . These two-level atoms interact with 
M bosonic modes, which have frequencies interact with coupling strengths {A Q }, and 

are described by the bosonic creation and annihilation operators {o^} an d { a a}- in terms 
of these quantities the full DH is given by 

N AI M N , 

H = u £ 4 l) + £ <*Afl« + E E 7= K + ««) (s? + s^) , (4) 

i=l a=l a=l i=l * 

where we have set h — 1. The origin of the factor 1/yN in the interaction is the fact that 
the original dipole coupling strength is proportional to 1/ VV, where V is the volume of the 
cavity. By writing p = N/V, where p is the density of the atoms in the cavity, this becomes 
■sj p/N and by subsuming the density into the coupling constants, {X a }, we obtain 1/y/N 
explicitly in the coupling. 

In Eq. (|3J) we have not made the usual rotating-wave approximation (RWA) under which 
one would neglect the counter-rotating terms a a s^} and a a s_ . We shall consider aspects of 
of the RWA in section EU 

We now specialise the Hamiltonian to consider a single mode bosonic field, and thus 
we drop the subscript a. The analysis of this Hamiltonian is further simplified by the 
introduction of collective atomic operators, 

N N 

j^E*?; ^ = E*±- ( 5 ) 

i=i i=i 

These operators obey the usual angular momentum commutation relations, 

[Jz,J±] = ±J±; [J+,J-] = 2J Z . (6) 

The Hilbert space of this algebra is spanned by the kets 
{\j,m); m — — j, — j + 1, . . . ,j — 1, j}, which are known as the Dicke states, and 
are eigenstates of J 2 and J z \ J z \j,m) = m\j,m) and J 2 |j, m) = j(j + l)\j,m). 
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and in terms o 
entanglement 



The raising and lowering operators act on these states in the following way: 
J±\j,rn) = a/j (j + 1) — m (m ± 1) \j,m ± 1). Note that j corresponds to Dicke's 
"co-operation number" which takes the values |, |, . . . , y for A" odd, and 0, 1, . . . , y for 
iV even. For example, with N = 2 atoms, j can take the values and 1. In terms of 
the s z values of the individual spins, the sector with j = 1 contains the triplet states 
ID) ^75 (I Tl) + I IT)) an d I TT)- The j = sector contains only the singlet state, 
"^72 (I TX) — I XT))- 111 general, the set of atomic configurations for iV > 2 is non-trivial j^jj], 
the individual atom configurations, the states are non-separable and contain 
4l| . In this work, we shall take j to have its maximal value, j = N/2, and 
once set, this value of j is constant, as the interaction in the DH does not mix j-sectors. 
Thus, the collection of N two-level systems is described as a single (N + l)-level system, 
which is viewed as a large pseudo-spin vector of length j = N/2. 

In terms of the collective operators, the single-mode DH may be written 

H = lu J z + ua ] a + —= (a j + a) ( J+ + J_) . (7) 
V 2j 

In the following, when we refer to the Dicke Hamiltonian we shall mean this single-mode 
Hamiltonian unless otherwise stated. The resonance condition is u = uq, and when plotting 
results we generally work on scaled resonance, such that u) — Uq — 1. 

Associated with the DH is a conserved parity n, such that [H, U] = 0, which is given by 

n = exp |ivrA>| ; N = a f a + J z + j, (8) 

where iV is the "excitation number" and counts the total number of excitation quanta in 
the system, n possesses two eigenvalues, ±1, depending on whether the number of quanta 
is even or odd, and correspondingly the Hilbert-space of the total system is split into two 
non-interacting sub-spaces. 

If we express the Hilbert-space of the total system in terms of the basis {\n) <g> \ j,m)}, 
where |n) are number states of the field, a^a\n) = n\n), and \j, m) are the Dicke states, the 
DH and the significance of the parity operator may be viewed in a simple lattice analogy. 
We construct a two-dimensional lattice, each point of which represents a basis vector and is 
labeled (n, m) . An example of this lattice with j = 1 is shown in Fig. ^ Note that the lattice 
is finite in the 'm' direction, but infinite in the 'n' direction, reflecting the dimensionality of 
the Hilbert-space. In this picture, we see that because the interaction conserves the parity 
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n 



FIG. 1: Schematic lattice representation of the states of the Dicke model for the example of j = 1. 
Shaded (unshaded) dots denote states of positive (negative) parity, with solid lines representing 
the couplings between the states. 

IT, states with an even total excitation number n + m+j interact only with other even states, 
and odd states interact only with odd states. This has the effect of dividing the total lattice 
into two inter-weaved sub-lattices, which correspond to the two different parity sectors. 
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III. THERMODYNAMIC LIMIT 



We begin by considering the DH in the thermodynamic limit, in which the number of 
atoms becomes infinite, N — > oo, and hence j — > oo. In this limit, the DH undergoes a 
QPT at a critical value of the atom-field coupling strength A c = ^ujuj /2, at which point 
the symmetry associated with the parity operator II of Eq. (JSj) is broken. To describe this 
QPT we shall derive two effective Hamiltonians, one to describe the system in the normal 
phase A < A c , and one to describe it in the broken-symmetry, super-radiant phase A > A c . 
It should be noted that the results derived below are exact in this limit, and this allows us 
to understand the nature of this system in a very detailed way. 

In this analysis we shall make extensive use of the Holstein-Primakoff representation of 
the angular momentum operators, which represents the operators in terms of a single bosonic 



mode in the following way 



29 



3, 



J+ = b ] ^/2j-tfb ■ J_ = y/2j-tfb b 

Jz = {tfb - j) , (9) 

where the introduced Bose operators obey [6, 6^] = 1. 

Making these substitutions into the DH of Eq. ((7j), we obtain the two-mode bosonic 
Hamiltonian 



H = lu (tfb - j) + ujcJa + A (a f + a) 1 - ^ + y 1 - ^bj . (10) 

In this representation the parity operator n becomes 

n = exp {in [a ] a + b ] b] } , (11) 
and the analogy with the standard parity operator of a two-dimensional harmonic operator 



is thus apparent 



42] 



A. Normal phase 

We derive an effective Hamiltonian for the system in the normal phase by simply neglect- 
ing terms in the full Hamiltonian of Eq. (fTU|) with j in the denominator. This approximates 
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the square-root in the Holstein-Primakoff mapping with unity, and we obtain the effective 
Hamiltonian given by 

H {1) = u tfb + uala + A (a f + a) (6 f + 6) - jcuo, (12) 

which is bi-linear in the bosonic operators and can thus be simply diagonalised. This is 
most easily facilitated by the introduction of position and momentum operators for the two 
bosonic modes, 

x = — = (a f + a) ; p x = iyf\ (a f - a) 

V 2to> 



y = -i=(&t + 6 ) ; J^i^f (tf-fc). (13) 

This representation will be particularly useful when we come to consider the wavefunctions 
of the system. Expressing Hamiltonian in terms of these operators we obtain 

if M = 1 { w V + p 2 + cj 2 y 2 + p 2 + 4Av^ xy - uj - lo} - ju , (14) 

which may be diagonalised by rotating the coordinate system in the following way 

x — qi cos 7^) + q 2 sin 7^; y = — q 1 sin 7^ + q 2 cos 7^^, (15) 

where the angle 7^ is given by 

tan(2 7 «H^. (16) 

On resonance, u = u , 7^ = 7r/4, so that x = (gi + q 2 ) /V% and y = (—q± + #2) /\/2- This 
rotation eliminates the xy interaction term in the Hamiltonian, which then assumes the form 
of two uncoupled oscillators, 

H {1) = \ {e ( ^q\ + V \ + efq 2 2 + p 2 2 - u - u ) - ju . (17) 

We now re-quantise with the introduction of two new bosonic modes defined by 

Qi = 7=== (4 + c i) ; pi = i 4\ {A ~ Cl 



2e W 



Q2 



= (4 + c 2 ) ; p 2 = i\f^- (4 - c 2 ) , (18) 



2e (1) 
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and arrive at the final diagonal form 

# (1) = e^c\ Cl + e^c\c 2 + \ ( e « + s W - u, - u ) - ju> . (19) 

The bosonic operators jci, c\, C2, <4 j, in terms of which is diagonal, are linear combina- 
tions of the original operators {a, a\ b, as detailed in appendix^! and describe collective 
atom-field excitations. The energies of the two independent oscillator modes are given 
by 

4 1)2 = \ + ± \l(co 2 -u 2 ) 2 + 16\ 2 couo^ . (20) 

Crucially, we see that the excitation energy e_ is real only when uo 2 + uo 2 > 
J(a;2-^ 2 )' + 16A 2 tutuo; or equivalently A < = A c . Thus we see that i^ 1 ' re- 

mains valid for A < A c , i.e. in the normal phase. In this phase, the ground-state energy is 
given by Eq = — jco> , which is O(j), whereas the excitation energies e± are 0(1). This 
means that scaling our energies with j, the excitation spectrum above the ground state 
becomes quasi-continuous in the j — > oo limit, that is to say that the excitation energies 
differ by an infinitesimal amount from Eq. 

It should be noted that H^ l > commutes with the parity operator II, and thus the eigen- 
states of H( l > have definite parity, with the ground state having positive parity. This can 
been seen from the fact that at A = 0, the ground-state is |0)|j, — j) in the original \n)\j, m) 
basis, which clearly has an even excitation number, n + m + j = 0. As the energy levels in 
the normal phase are non-degenerate, the continuity of the ground state with increasing A 
ensures that it always has positive parity in this phase. 



B. Super-radiant phase 

In order to describe the system above the phase transition, we must incorporate the fact 
that both the field and the atomic ensemble acquire macroscopic occupations. To do this, 
we start with the Holstein-Primakoff transformed Hamiltonian of Eq. (fTTH) and displace the 
bosonic modes in either of the following ways 

a f -> c f + v^; b ] -> d ] - <J@. (21) 

or 

a f -> c f - v^; b ] -f S + y/p. (22) 
12 



Crucially, we assume that the as yet undetermined parameters a and j3 are of the O(j), 
equivalent to assuming that both modes acquire non-zero, macroscopic mean-fields above 
A c . In the following we shall just consider the displacements given by Eq. (|21|). as the 
calculation with the other choice is identical but for a few changes of sign. 
Making these displacements, the Hamiltonian of Eq. (|1U|) becomes 

H = u |rf f d - (d j + d) + (3 - j\ + uo {c j c + y/a (c f + c) + a} 

+ A^(c t + C + 2 V ^) (fy/i+y/ld-2y/Py/t), (23) 



where for brevity we have written = y 1 — ~ an d k = 2j — /3. Taking the 

thermodynamic limit by expanding the square-root y/£ and then setting terms with overall 
powers of j in the denominator to zero, we obtain 



H i2) = + \u + y W ^ K d ~ I 2A J ^ " ^ \ (c f + c) 



+ <!c (/3-j)+^-^^(l + 4A;)J>. (24) 

We now eliminate the terms in the that are linear in the bosonic operators by choosing 
the displacements a and (3 so that 

2\J^-u^=0, (25) 

and 

— U-P)-uo)y/0 = O. (26) 

The a//3 = \/a = solution of these equations recovers the normal phase Hamiltonian 
The non-trivial solution gives 



« = ^|(i-/i 2 ) ; v^VTTw), (27) 
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where we have defined 



H = 



LUUJq A 2 



AX 2 A 2 ' 

With these determinations, the effective Hamiltonian of Eq. (|24|) becomes 



(28) 



^_ ctc+ | (1 + , )dt , + Ml_^ (d t + d): 



_|_ (C , +C)( „ + , ) _ J {^ + ^ 



(29) 



To facilitate the diagonalisation of this bi-linear Hamiltonian we move to a position- 
momentum representation defined by 



X = _=( C t + C ); p x =,y|( c t_ c ) 
v 2cj 



x 2w 



+ tZ) ; Py 



(rft _ d) 



(30) 



where a) = ^ (1 + (/,). Note that this is not the same representation as defined in Eq. ()13|). 
The diagonalisation then proceeds similarly to before, involving a rotation in the X-Y plane 
to the new coordinates 



X = Q 1 cos 7 (2) + Q 2 sin 7 (2) 
Y = -g lS in7 (2) + Q 2 cos7 {2) 



(31) 



with the angle is given by 



tan(2 7 ( 2 ))= ^°f 2 . 

UJq — /Z 2 UT 



(32) 



Subsequent requantisation in terms of two new modes, e± , corresponding to the rotated, 
decoupled oscillators gives us the diagonal form 



1 I ( 2 ) , ( 2 ) W 0/1 \ 

+ 2T V + - 1 + 

with the oscillator energies being given by 



u (1- n) 



(33) 



2ei 2)2 = ^ + ^± 



+ 4w 2 w 2 



(34) 
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The Bogoliubov transformations that induce this diagonalisation are given in appendix 
1X1 The excitation energy and hence H^ 2 \ remains real provided that ^% + u 2 > 



A 4 



,2 

^ -LU 2 



2 

+ 4lj 2 Uq, or equivalently A > = A c . Thus we see that describes 

the system in the super-radiant phase, A > A c , in which the scaled ground-state energy is 
given by4 2) /j = -{^ + ^}. 

If we choose the signs of the operator displacements as per Eq. we obtain exactly 

the same values of a and (3, and an effective Hamiltonian identical in form with Eq. (}33|) . 
This clearly has the same spectrum and therefore, each and every level of the total spectrum 
is doubly degenerate above the phase transition. What has occurred is that the symmetry 
of the ground state, defined by the operator II, has become spontaneously broken at A c . The 
Hamiltonian H^ 2 \ for either choice of displacement, does not commute with II, and thus its 
eigenf unctions do not possess good parity symmetry. 

Although the global symmetry II becomes broken at the phase transition, two new local 
symmetries appear, corresponding to to the operator 

II (2) = exp {m [c ] c + Sd] } , (35) 

for both the two different choices of mean-field displacements. This operator commutes with 
the appropriate super-radiant Hamiltonian, [if( 2 ) , IB 2 '] = 0. 

C. Phase Transition 

Having derived the two effective Hamiltonians which describe the system for all A in 
the j —* oo limit, we now describe the systems properties in each of its two phases. The 
fundamental excitations of the system are given by the energies e±, which describe collective 
modes, similar to polariton modes in solid-state physics [43 1. The behaviour of these energies 
as a function of coupling strength is displayed in Fig. |2j where we have labeled the two 
branches as "atomic" and "photonic", according to the nature of the excitation at zero 
coupling. From this figure we see that as the coupling approaches the critical value A c , the 
excitation energy of the photonic mode vanishes, £_ — * 0, as A — > A c , demonstrating the 



existence of the QPT. In contrast, e + tends towards a value of a/^o + u;2 as A — > A c from 
either direction. In the asymptotic limit of A —>■ 00, s_ — > uj (returning to its A = value) 
whereas e + — > AX 2 /uj. The critical exponents of this QPT are manifested in the behaviour 
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FIG. 2: The excitation energies of the Dicke Hamiltonian in the thermodynamic limit as a function 
of coupling A. The Hamiltonian is resonant, u = uq = 1, and the vanishing of e_ at A = A c = 0.5 
signals the occurrence of the QPT. 

of the excitation energies Q|. As A — > A c from either direction, the energy e_ can be shown 
to vanish as 

£ ( A _» Ac) „ / 3 ^ IA C - \\ 1 ' 2 . (36) 

v ; y 16A4+0J 4 ' 1 v ; 

The vanishing of at A c reveals this to be a second-order phase transition. We define the 
characteristic length scale in the system in terms of this energy as 

I- = 1/y/eZ. (37) 

From Eq. ()36)1 this length diverges as | A — A c | ^ with the exponent v — 2. We then write 
that E- vanishes as |A — \ c \ zv ', with the dynamical critical exponent being given by z = 2. 
At the phase transition point, we have 

(A c ) = # (2) (A c ) = Vlu 2 + w 2 4c2 + l - [Vuj 2 + to 2 - u - o; ) - ju , (38) 

from which we see that at A c the system becomes effectively one-dimensional. 

The ground-state energy of the system Eq is shown in Fig. El and the analytic form 
expression is given in Table HI Note that we scale all quantities by j, which means that 
the plotted Eq/ j is equal to 2Eq/N, twice the energy per atom. We also plot the second 
derivative of the ground-state energy with respect to A, which possesses a discontinuity at 
A c , clearly locating the phase transition. 
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FIG. 3: The scaled ground-state energy Eq/j and its second derivative j~ 1 d 2 EQ/d\ 2 as a func- 
tion of coupling A. Solid lines denote results in the thermodynamic limit, whereas dashed lines 
correspond to the results for various finite values of j = 5, 1, |, 3, 5. The Hamiltonian is resonant: 
lj = u)q = 1, A c = 0.5. 





A < A c 


A > A c 


E G /j 
(J z )/j 


-u 
-1 



oj X' 2 uj 

-a c 2 M 2 

2 (A 4 - A 4 ) / („A) 2 



TABLE I: The ground-state energy, atomic inversion and mean photon number of the Dicke 
Hamiltonian in the thermodynamic limit. 

In Fig. 0| we plot the atomic inversion (J z )/j and the mean photon number (n a )/j = 
(a'a)/j. This figure clearly illustrates the nature of the phase transition, - in the normal 
phase, the system is only microscpically excited, whereas above A c both the field and atomic 
ensemble acquire macroscopic excitations. We may write the values of the atomic inversion 
and the mean photon number above A c in the following fashion: 

(J z )/j = 1 - 13/ j , {a} a)/ 3 = a/ j ; A > A c . (39) 

Thus making clear the physical meaning of the displacement parameters a and (3 of Eqns. 

flZZj). 
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FIG. 4: The scaled atomic inversion and mean photon number of the Dicke Hamiltonian as a 
function of coupling A. Solid lines denote results in the thermodynamic limit, whereas dashed lines 
correspond to the results for various finite values of j = ^, 1, |, 3, 5. The Hamiltonian is resonant: 
lj = loq = 1, A c = 0.5. 

D. Ground-state Wavefunction 

We now consider the ground-state wavefunctions of the system above and below the 
phase transition. After diagonalisation, the two effective Hamiltonians are both of the 
form of a pair of uncoupled harmonic oscillators. Thus, in the representation in which the 
Hamiltonians are diagonal, their wavefunctions will simply be the product of the appropriate 
harmonic oscillator eigenf unctions. Here we seek to express these wavefunctions in terms 
of the two-dimensional x-y representation of Eq. (|13|) - which corresponds to the original 
atomic and field degrees of freedom. 

We have already noted that in the Holstein-Primakoff representation the parity operator 
las the form n = exp [in + From our knowledge of the harmonic oscillator 

. we know that the action of n in the x-y representation is to perform the coordinate 
inversions, x — > —x and y — > — y, with p x and p y remaining unaffected. Thus the operation 
of n is equivalent to a rotation of it about the coordinate origin and, in the normal phase 
where n is a good quantum number, the wavefunctions will be seen to be invariant under 
this rotation. 

The ground-state wavefunction of a single harmonic oscillator in terms of its coordinate 
q is a Gaussian with width determined by the energy of the oscillator. Correspondingly we 
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FIG. 5: The ground-state wavefunction of the low-coupling Hamiltonian in the x-y 
position-momentum representation for couplings A = 0, 0.3, 0.49, 0.4999999. The Hamiltonian is 
resonant: u = ujq = 1, A c = 0.5. 

define the normalised Gaussian functions 

/ (1,2)\V4 f (1,2) 1 

M exp|-^-g 2 L (40) 

where e^' 2 ^ are the excitation energies encountered earlier. 

In the normal phase, the effective Hamiltonian is diagonal in the q±-q2 representation 
of Eq. and its ground-state wavefunction Wq in this representation is therefore 

*S ) (ft,<6)=GL 1) (ff 1 )Gi 1) (< & ). (41) 

Moving to the x-y representation, we have 

^ (x, y) = G W (x cos 7 (1) - y sin 7 {1) ) G$ (x sin 7 (1) + y cos 7 {1) ) , (42) 

and this wavefunction is plotted for various couplings in Fig. At A = the wavefunction is 
the product of orthogonal Gaussians of equal width (on resonance). As coupling increases, 
the wavepacket becomes stretched in a direction determined by the angle 7^ 1 - ) , which on 
resonance is simply equal to 7r/4. This stretching increases up to A c , where the wavefunction 
diverges. We thus see the significance of the length /_ introduced earlier - it is the extent 
of the wavefunction in the direction of this stretching. Correspondingly, l + is the extent of 
the wavefunction in the orthogonal direction. 

In the super-radiant phase, the ground state is degenerate. We shall initially consider 
the ground-state wavefunction of the effective Hamiltonian with displacements chosen 
in Eq. (j'21j) . This is diagonal in the Q1-Q2 representation of Eq. (jUlj) . and therefore its 
ground-state wavefunction is 

^g ) (Qi,g 2 ) = G (2) (g 1 )Gf(Q 2 ). (43) 
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FIG. 6: The ground-state wavefunction of the high-coupling Hamiltonian in the x'-y' 
position-momentum representation for couplings A = 0.500001,0.51,0.6,1.0. The Hamiltonian is 
resonant: u = ojq = 1, A c = 0.5. 

Using Eqns. (pHIj) . (}2*T]) . (jSOjl. and (|3Tj) we may write this in the original x-y representation 
as 

(x,y) = G {2) (Jx - a/M^) cos 7 (2) - (y + VWJ^) sin 7 (2) ) 

xGf ((a: - sin 7 ^ + (y + y/2~P/^ cos 7 (2) ) . (44) 

This expression contains displacements involving the macroscopic quantities a and /3, and 
so we define the new coordinates x' and y' to remove them: 

x' = x - A x ; y' = y + A y , (45) 

with 

A x = y/2a/u; A y = y/2/3/u , (46) 

which are both proportional to \fj. The relationship between the coordinate system X-Y 
of Eq. (|3Ujl and x'-y' is very simple, namely x' — X and y' = \fuj/ujQ Y . The coordinate 
system x'-y' is useful because although X-Y is the diagonal representation for the super- 
radiant phase, the definition of these coordinates depends upon Co and hence upon A, which 
distorts the picture. In terms of these coordinates the wavefunction becomes 

*g> (x',y') = (co /uj) 1/4 (Vcos 7 (2) - y^/tiy'sm^ 

xGf (Vsin 7 (2) + v 4V^y , cos 7 {2) ) . (47) 

Figure El shows W G (x',y') for four different couplings. Just above the phase transition 
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FIG. 7: Parametric plot of the (scaled) displacements (A x ,— Ay) and — (A x ,A y ) as A is varied 
between 0.5 and 10. The Hamiltonian is resonant; u> = loq = 1. 

the wavefunction is in a highly deformed state, characterised by the divergent Z_. As the 
coupling increases further above A c , the wavefunction relaxes back to a well localised state. 

(2) 

When considered in the original x-y representation, the wavefunction ty G (x', y') pictured 
in Fig. El is centered about the point (+A X ,— A y ), which lies in the lower-right quadrant 
of the x-y plane. The complementary wavefunction, identical in shape with this one but 
determined by the displacements is centered at (— A x , +A y ), in the upper-left quadrant. 
The positions of these two centers as parametric functions of coupling are shown in Fig. [71 
These two wavefunctions, corresponding to the two choices of displacement, are separated 
from the origin by an amount proportional to \fj. It is thus clear that neither of these 
wavefunctions is symmetric under rotation of it about the origin of the x-y coordinate 
system, demonstrating once more that the II symmetry has been broken. There is, however, 
symmetry with respect to a rotation of it about the origin of each x' — y' coordinate system, 
which corresponds to the existence of the local symmetries associated with of Eq. (|35j). 

It is interesting to consider the behaviour of the ground-state wavefunction as A —>■ oo. In 
this limit — > oj, e^} 4X 2 /uj and the mixing angle of the two modes 7*- 2 - 1 tends to zero, 
meaning that the modes decouple. The Bogoliubov transformations of the modes become 



(3d f + d) 



e 2 



+ 3d) 



(48) 



2y/2 v ' 2V2 

illustrating the decoupling. Note that the ex simply reverts to the c mode, whereas the e 2 

mode tends towards a linear combination of the annihilation and creation operators. In this 
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limit, the wavefunction becomes 



— exp (^-cjoZ/ - -x ) , (49) 



which is independent of A. 



E. Squeezing 

A bosonic field may said to be squeezed if the uncertainty in either of its quadratures, (x 
or p x ) is less than the uncertainty in a coherent state Q]. A coherent state is a minimum 
uncertainty state with (Ax) 2 (Ap x ) 2 = 1/4 and with the uncertainty apportioned evenly 
between the two quadratures. Therefore, the field is squeezed whenever (Ax) 2 or (Ap^) 2 



has a value lower than 1/2 |45[ . 

We define the two quadrature variances of the original field mode a by (Ax) 2 = (x 2 ) — (x) 2 
and (Ap x ) 2 = (p 2 ,) — (p x ) 2 , which may be seen to be equal to 

(Ax) 2 = i- {l + (a t2 ) + (a 2 ) + 2(a t a) + ((a f ) + (a)) 2 } , 
(Ap x f = | {l - (at 2 ) - (a 2 ) + 2<ato> + ((at) + (a)) 2 } . (50) 

As we have introduced a bosonic algebra for the atomic collection, we now introduce an 
analogous definition for squeezing in the atoms, and say that in terms of the variances, 

(Ay) 2 = ^- {l + (&t 2 ) + (b 2 ) + 2 <&t 6) + ( (6 t) + (6) ) 2 | ( 
(A Vy f = y {l - (& f2 ) - (b 2 ) + 2(6t6> + «6t) + ( & )) 2 } . (51) 
the atoms are squeezed if either (Ay) 2 or (Ap y ) 2 is less than 1/2. The squeezing of atomic 



ensembles is usually defined in terms of the collective operators |46|,|47|. Because the angular 
momentum operators obey the commutation relation, [J + , J_] = 2J Z , the uncertainty rela- 
tion, (A J,;) 2 (AJ y ) 2 > \\{J Z }\ 2 holds for any state. By substituting in the Holstein-Primakoff 
forms into this expression and taking the thermodynamic limit, we see that this relation re- 
duces to (Ay) 2 (Apy) 2 > 1/4, demonstrating the equivalence in the thermodynamic limit of 
our definition in terms of y and p y and the usual one. 
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1,5 

FIG. 8: The squeezing variances of the ground state of the DH in the thermodynamic limit. The 
Hamiltonian is resonant: to = loq = 1, A c = 0.5. Note that on resonance, (Ax) 2 and (Ay) 2 are 
coincident for A < A c , and the same for (Ap x ) 2 and (Ap y ) 2 

In the normal phase the expressions for the variances are evaluated by simply making 
the appropriate substitutions from Appendix El and taking their ground-state expectation 
value. In the super-radiant phase, it can be shown that the variances of the original field and 
atomic modes of Eqs. (joTi|) and (JoTj) can be expressed in terms of the displaced coordinates 
as follows 

(Ax) 2 = (Ax') 2 = (AX) 2 ; (Ap x f = (Ap' x f = (AP X ) 2 , 

(Ay) 2 = {Ay') 2 = ^JI7 (AY) 2 ; (Ap y ) 2 = (Ap' y ) 2 = ^J^(APy) 2 . (52) 

This results from the fact that the squeezing variances do not depend upon the displacements 
of the field modes, and thus evaluating the super-radiant variances is as simple as in the 
normal phase. 

The analytic values of these variances in the ground state are shown in Appendix FBI and 
are plotted as functions of coupling in Fig. |HJ In the normal phase, as A approaches A c there 
is a sharp increase, and eventually a divergence, in (Ax) 2 and (Ay) 2 . This is accompanied 
by a slight squeezing of the momentum variances. In the super-radiant phase the initially 
divergent values of (Ax) 2 and (Ay) 2 reduce rapidly with increasing coupling. The behaviour 
of these variances reflects the nature of the wavefunctions plotted in Figs. EJandlHl Notice 
that as A — > oo, (Ax) 2 and (Ap x ) 2 return to their A = values, whereas (Ay) 2 and (Ap y ) 2 
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become squeezed and anti-squeezed respectively. This is in agreement with the results of Eq. 
([48)1 . which show that the e± mode becomes identical to the c-mode, which is unsqueezed, 
whereas the 62 mode reverts to a linear superposition of d and eft operators , w hich is a 
specific example of the Bogoliubov transformation producing a squeezed state 48^ . 
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IV. THE ONSET OF CHAOS 



As we have just demonstrated, the DH is exactly integrable in the thermodynamic limit. 
However, for finite j, this is not the case and the possibility of quantum chaos remains. The 
signature of quantum chaos that we use to investigate this possibility is the character of the 
energy spectrum as quantified by the nearest-neighbour level distribution P(S). Bohigas 
et al.p first conjectured that the study of spectral quantities such as P(S), and their 
comparison with the results from random matrix theory should give an indicator of quantum 
chaos. This may be understood by the following argument. Classically integrable systems 
have high degrees of symmetry and hence their quantum counterparts have many conserved 
quantum numbers. This permits level-crossings to occur in the spectrum, leading to a P (S) 
with a maximum at small level-spacing, S —>■ 0, with a P (S) given by the Poissonian 
distribution Pp (S) = exp (— S). We shall call quantum spectra with Poissonian statistics 
"quasi-integrable" . Conversely, classically chaotic systems have no such integrals of motion 
and we thus expect their quantum energy spectra to be highly correlated and absent of 
crossings, leading to P (S) — > as S —>■ 0. Although the precise form of the P (S) for chaotic 
systems depends on the symmetries of the model, we shall find that only the Wigner-Dyson 
distribution, Pw (S) = irS/2exp (—ttS 2 /4), is of relevance here 

Despite its popularity, it should be pointed out that the correspondence between the 
P (S) distribution and the intestability or otherwise of the classical system is not absolute, 

nn 

and exceptions do exist [50, 51]. Despite this, the P (S) does provide a convenient and useful 
indicator of quantum chaos, and the conjecture does hold true in a countless examples. In 
the case in hand, this signature turns out to be very accurate, as will be evinced when we 
compare the P (S) results with those of our semi-classical model. 

A. Numerical diagonalisation 

Exact solutions for the DH at finite j do not exist, except in the very special case of 
j = 1/2 where isolated exact ("Juddian") solutions may be found 0, [^J . Consequently we 
employ numerical diagonalisation to investigate the system. To perform these diagonalisa- 
tions we use the basis {\n) ® \j, m)}, where \n) are number states of the field, and \j, m) are 
the Dicke states. In performing the diagonalisation, we truncate the bosonic Hilbert space 
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but always maintain the full Hilbert space of the pseudo-spin. The size of the matrices re- 
quiring diagonalisation is reduced by restricting ourselves to a single parity subspace, which 
is achieved by only considering states with n + m + j even or odd for positive and negative 
parity respectively. With j finite, IT is a good quantum number, independent of coupling, 
and the ground state always has positive parity. 

The results obtained via this diagonalisation for the ground-state energy and its second 
derivative for a sequence of finite j values have been plotted alongside the j ; — > oo results 
in Fig. |3J whilst the corresponding atomic inversions and mean photon numbers are plot- 
ted in Fig. HI These figures demonstrate how rapidly the finite j results approach their 
thermodynamic limits as j is increased. 



B. Level statistics 



Having numerically obtained energy spectra of the DH, we can construct the nearest- 
neighbour level-spacing distribution P(S). This is formed from a large number of levels 
from the spectrum, which we initially unfold to rid of secular variation [49]. We then 
calculate the level spacings, 

S n = E n+ \ — E n , (53) 

where {E n ; n — 0, 1, . . .} is the set of eigen-energies of the DH with positive parity, and 
construct their distribution function P (S). Finally, we normalise the results for comparison 
with the universal ensembles of Random Matrix Theory j3| . 

Figure El shows the P (S) distributions obtained for the DH at various values of A and j. 
At low j (j < 3) the P (S) clearly do not correspond to any of the universal ensembles, but 
rather to non-generic distributions consisting of several isolated peaks. This is most obvious 
in the j = 1/2 case (not shown here), which is known as the Rabi Hamiltonian (RH) [54]. 
The RH has a spectrum that is of "picket-fence" character j^, which is characteristic 
of genuinely^ integrable models such as one-dimensional systems and systems of harmonic 



oscillators P). The RH is unusual and must be treated as a special case because, although it 
has never been shown to be integrable, isolated exact solutions do exist |52u53|. Moreover 
the model is separable and may be reduced to a single degree of freedom |32j. 

Returning to the P (S) distributions, we see that at low couplings A < A c , (for example, 
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FIG. 9: Plots of nearest-neighbour distributions P(S) for the Dicke Hamiltonian, for different 
couplings A and pseudo-spin j. Also plotted are the universal Poissonian (dots) to Wigner (dashes) 
distributions. The Hamiltonian is resonant: u = Uq = 1, A c = 0.5. 

A = 0.2 in Fig. @), as we increase j, the P (S) loses its non-generic features and approaches 
ever closer the Poissonian distribution, Pp (S). At and above the critical coupling (A = 0.5 
and 0.8 in Fig. |HJ) the spectrum is seen to converge onto the Wigner-Dyson distribution 
Pw (S) as j is increased. 

The nature of the change in the P (S) distribution may be characterised by the quantity 

jf [P (S) - P W (S)} dS 



n 



f So [P P (S) - P W ds 



(54) 



where Sq = 0.472913 . the value of S at which the two generic distributions Pp (S) and 
Pw (S) first intersect rj measures the degree of similarity of the calculated P (S) to the 
Wigner surmise Pw (S), and is normalised such that if P (S) = Pw {S) then r] = 0, and if 
P (S) = Pp (S) then rj — 1. The behaviour of rj as a function of coupling for j = 5 and 
j = 20 is shown in Fig. Considering the j = 20 case first we see that the spectrum 
is strongly Poissonian at low couplings and that at A is increased towards A c , it becomes 
more Wigner-Dyson like. This proceeds until we reach A c , about which the spectrum is 
remarkably well described by Pw (S). Note that for A < A c the value of rj drops steadily 
with coupling, whereas above A c it maintains an approximately constant value close to zero. 
For the j = 5 case, a similar transition is observed, but it is not as pronounced and the 
agreement with the universal distributions is not as good as in the higher j case. 
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FIG. 10: The modulus of 77, Eq. (|54|). plotted as a function of coupling for systems of j = 5 and 
j = 20. A value of r\ = 1 indicates Poissonian statistics and r] = corresponds to Wigner-Dyson. 
The system is on scaled resonance {u = Uq = 1). 

Thus, for sufficiently high j, we see a significant change in P (S) as A is increased from 
zero through the critical value A c . Below A c , there is a significant amount of level-crossing, 
which decreases as we approach A c . Above A c there is practically none, to within statistical 
error. Thus we conclude that the precursors of the QPT in this model lead to a cross-over 
from quasi- integrable to quantum chaotic behaviour at A ~ A c for sufficiently high j. 




C. Regularity at low energy 

A further transition between integrable and chaotic behaviour is observed in the sequence 
of level spacings S n as the coupling is increased from A c to oo. In the A — > oo limit, the DH 
is integrable for arbitrary j and equivalent to 

-f^A^oo = uja^a + 2—^= (a f + a) J x . (55) 

The eigenstates of H x ^oo ar e obviously eigenstates of J x , and thus 

Hx^oo = ua)a + 2m— — (at + a) , (56) 
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FIG. 11: Nearest-neighbour spacing S n = E n+ \ — E n vs. eigenvalue number n plot for j = 5 with 
A = 4. Horizontal crosses: results for the integrable A — > oo Hamiltonian. Inset: j = 5 results with 
A = 2 and A = 3. The Hamiltonian is resonant: u = loq = 1, A c = 0.5 



where m = —j, . . . ,j is the eigenvalue of J x . This bosonic Hamiltonian is diagonalised by 
the displacement a — > a — 2mA/^/2j, giving the eigenvalues to be 

J ^3 

where k — 0, 1, 2, . . .. The energy levels with +m and — m are degenerate. 



(57) 



As A is increased from A c to approach this A — > oo limit with j fixed, the spectrum reverts 
from Wigner-like to integrable. However, it does not follow the usual transition sequence 
of Wigner distribution gradually changing into a Poissonian one, as one might expect, but 
rather through a sequence illustrated by Fig. El For couplings sufficiently higher than A c , 
the spectrum becomes very regular at low energy, where it approximates the A — ► oo of 
results very closely. Outside the regular region the spectrum is well described by the Wigner 
surmise, and the energy-scale over which the change between the two regimes occurs is seen 
to be surprisingly narrow. As coupling is increased, the size of the low-energy integrable 
window increases, until it eventually engulfs the whole spectrum as A — ► oo. This division 
of the spectrum into regions is close to Percival's conception of how regular and irregular 
behaviour would manifest itself in quantum systems 
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D. Wavefunctions for finite j 

We now consider the wavefunctions of the DH at finite j. To do this, we shall use the 
position-momentum representation of Eq. (J13j) used earlier in discussing the wavefunctions 
in the thermodynamic limit. We begin with the eigenfunctions obtained from numerical 
diagonalisation, which are of the form 

|*nm> = £ E C 2»U,™>, (58) 
n=Q m=—j 

where n c is the maximum boson number in the artificially truncated Fock space, and Cnm 
are coefficients. The position representatives of the number states of the field |n) are simply 
the usual Harmonic oscillator eigenfunctions 

( X \n) = ^^e-^ 2 H n (V^x) (59) 

where H n is the nth Hermite polynomial. For the angular momentum part of the basis vector, 
we recall that under the Holstein-Primakoff mapping J z — > b^b — j, and thus the Dicke states 
are eigenstates of tfb with eigenvalue (j + m): b^b\j,m) = (j + m) \ j,m), —j < in < j . 
Consequently, we may represent the Dicke states in the same way as the Fock states above, 
allowing us to write the total wavefunction in the two-dimensional position representation 
as 

{x,y) = V^VK-W) £ V c ^if^r^ - (60) 

nm\ ,y/ ^ Z— ✓ nm 2( n +-?+ m ) Til ( 1 '+ Wl) ! 

n=0 m=~j " ^ J ' 

This is a very productive representation in which to study this Hamiltonian. It does however, 
suffer from the drawback that whereas the set of oscillator eigenfunctions Eq. (}59|) forms an 
orthonormal set in the x direction, this is not the case in the y direction as we only keep 
up to the (2j)th oscillator eigenfunction in this direction. This means, for example, that 
we could not go from an arbitrary wavefunction in the y direction to a description in terms 
of the Dicke states, because we do not have a complete set of functions in this direction. 
Specifically, the significant width of the wavefunction is limited in the y direction by the 
maximum significant extent of the highest Hermite polynomial i^j- However, if we know 
the value of j and only consider wavefunctions that are describable in terms of these then 
the representation is unique. 



30 



• 









-4.5 4.5 -4.5 4.5 



FIG. 12: The modulus of the ground-state wavefunction ip(x,y) of the Dicke Hamiltonian in 
the abstract x-y representation for finite j = 5, at couplings of A/A c = 0.2, 0.5, 0.6, 0.7. Black 
corresponds to Max|^| and white corresponds to zero. The Hamiltonian is resonant u> = ujq = 1; 
A c = 0.5. 

Figure [T21 shows the ground-state wavefunction of DH with j = 5, for a series of increasing 
couplings. Note that the wavefunction is always invariant under a rotation of tt about the 
origin as demanded by the II symmetry. This wavefunction starts as a single lobe centred at 
the origin for low coupling. As the coupling increases, the two modes start mixing, leading 
to a stretching of the single-peaked wavefunction, which then splits into two as coupling is 
increased through a coupling approximately equal to A c . With further increases in coupling 
the two lobes move away from each other in their respective quadrants of the x — y plane. 

The key observation regarding the two lobes formed above A c is that, provided j and A 
are both sufficiently large, their displacement from each other is proportional to y/j, and 
that this is a macroscopic quantity. The excited states exhibit similar behaviour, having an 
extent proportional to \fj above the phase transition. 

Therefore, around the critical coupling A ~ A c , the wavefunctions of the system become 
delocalised, and the extent of this delocalisation is proportional to \fj. As this is a macro- 
scopic quantity, we see that above A c , the system at finite j develops macroscopic coherence 
in its wavefunctions. The most striking example of this is the ground state, where the two 
macroscopically different lobes are reminiscent of the two states of a Schrodinger's cat. 

The delocalisation and accompanying macroscopic coherence are rather general features 
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of the onset of chaos, and are natural consequences of the exponential divergence of trajec- 
tories in a classically chaotic system 3]. If we consider a small volume of initial conditions 
in the classical phase space (a well-localised initial wave-packet), and let the system evolve 
chaotically, this initial volume rapidly becomes smeared out over the entire phase-space 
accessible to it. This is reflected in the quantum system by the delocalisation of the wave- 
functions. That such systems are macroscopic coherent may be seen from the observation 
that under Hamiltonian dynamics, the volume of the initial "wave-packet" remains constant 
in time (Liouville's theorem). This means that the exponential divergence in some direction 
leads to the exponential contraction in others. This contraction will continue until the size 
of the packet becomes of the order of h and quantum effects come into play. If we imagine 
that the wave-packet becomes narrow in the direction of momentum p, then the uncertainty 
Ap becomes very small. In order that the Heisenberg uncertainty relation holds, the uncer- 
tainty in the corresponding coordinate Ax must become very large, and this leads to the 
emergence of macroscopic coherence in the system. 

This effect is what was observed in the variances calculated earlier in connection with 
squeezing in the thermodynamic limit. As A — * A c from below, the variances (Ax) 2 and {Ay) 2 
diverged, with (Ap x ) 2 and (Ap y ) 2 remaining near their quantum limit of 1/2. The behaviour 
of these variances then reflects the onset of quantum chaos and the macroscopic coherence 
of the wavefunctions. A vital difference between the j — ► oo and the finite j results thus 
emerges in the super-radiant phase. In the thermodynamic limit, the variances (Ax) 2 and 
(Ay) 2 reduce as A is increased from A c , indicating that the wavefunctions become localised 
and lose this macroscopic coherence. Contrast this with the finite j, where sufficiently above 
A c the wavefunction is always delocalised and the variances are O (\/J). This is because, 
whereas at finite j we obey II symmetry and thus have both of the lobes of the wavefunction, 
in the thermodynamic limit, we consider each lobe separately under the broken symmetry. 
The individual lobes are themselves localised and this is where the discrepancy comes from. 
This is, we believe, the reason why although the spectrum is of the Wigner-Dyson type for 
large j, the spectrum in the j — > oo limit is integrable, as in this limit the wavefunctions 
possess no delocalisation and no macroscopic coherence. 

This picture also provides us with an explanation of why the P (S) for very small j are of 
the non-generic one-dimensional type. As the extent of the wavefunction in the y-direction 
is effectively constrained by the number of harmonic eigenf unctions in that direction, which 
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is determined by j, having a small j prevents full delocalisation in this direction, inhibiting 
the chaoticity of the quantum system. 
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V. THE SEMI-CLASSICAL MODEL 



As noted in the introduction, there have been many different semi-classical models derived 



from the DH 



37, 
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39]. That there have been so many different approaches is a 
reflection of the fact that the quantum mechanical spin possesses no direct classical analogue. 
Nevertheless, semi-classical models can be constructed, and in the following we shall propose 
a new approach. Before this, let us briefly examine some of the previous work. 

A widely discussed approach is that of a Hartree-Fock type approximation in which one 



derives the Heisenberg equations of motion 
these equations by their expectation values 



'or the system and replaces the operators in 



321 ]. These are treated as classical variables 



and a set of non-linear equations of motion are obtained for them, which show classical 
chaos for certain parameter ranges |37J. Despite this, the above approach is not completely 
satisfactory as the motion only depends j in a trivial way. Furuya at al. have studied a 
classical model similar to the one we propose below [3^. They derived their semi-classical 
Hamiltonian by evaluating the expectation value of the DH in a state composed of a product 
of photonic and atomic coherent states, and this system was also shown to exhibit chaos. 
Despite the similarity of their model to ours, they did not discuss the role of the phase 
transition in determining the chaoticity of the model, which is a key feature of our model. 
We start with the DH in the bosonic form of Eq. tjlOj): 



b ] J 1 - — + 

By using the inverse of the relations in Eq. (|13|) . namely 




(61) 



2 + ^ 



we may write this Hamiltonian in the position-momentum representation, 



UJ 

2 " 



— Po 



(62) 



H 



-j^o + - (u 2 x 2 + p 2 x -u + u 2 y 2 +p 2 y - u ) 



+ \y/iUUJ X 

where we have written 



V ~ —Py VT 

UJq 



n 



'1 



v[y + —p y 



(63) 



V = {^oV 2 +P 2 y~ Vo) I ( 4 J^o) 
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(64) 



We now move very naturally from this quantum-mechanical Hamiltonian to a semi-classical 
one by setting the position-momentum commutators to zero, i.e. = 0, [y,p y ] = 0. 

This causes the interaction term to become real and in terms of classical variables we have 

H sc = -JUo + \ (w V + pi ~ LQ + Ujlv 2 +P 2 y - Vq) 



, u%y 2 + p 2 -uj 

+ 2\^U^xy\ll — J . (65) 

Unusually, this Hamiltonian contains an intrinsic constraint, which is determined by the 
requirement that the square-root must remain real for the system to remain Hamiltonian. 
This means that the inequality 

^4^KV + ^-^o)<l (66) 

is satisfied for all time. 



A. Classical Phase Transition 



The Hamiltonian H sc undergoes a spontaneous symmetry-breaking phase transition that 
is directly analogous to the QPT of the quantum model. The exact correspondence between 
the classical and quantum Hamiltonians in the thermodynamic limit is because in this limit 
the system is exactly described with a mean-field theory as used earlier, and the use of 
classical variables as we have done here is equivalent to a mean-field theory. Consequently, 
we are able to derive classical effective Hamiltonians exactly as we did in the quantum case. 
The effective Hamiltonian for the normal phase is derived by simply letting j 1 — > oo (i.e. 
rj — > 0) in the Hamiltonian of Eq. (J65|) . This gives us 

H k = j + p x + u oV 2 +P 2 y + 4Av^) xy-uj Q -uj} - ju , (67) 

which is identical to Eq. from the quantum analysis, and may be diagonalised with 

the same rotation. The equilibrium position of Hamiltonian is the origin: x = y = 

Px =Py = 0. 

An effective Hamiltonian for the super-radiant phase is derived in the same way as in the 
quantum case, by displacing the co-ordinates as in Eq. so that x — > x'±A x ,y — > y'^fA y , 
where the displacements are the same as before: A x = ^J2a/uj and A y = ^J2(3/ujq. Making 
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these displacements and then taking the thermodynamic limit results in a Hamiltonian Hsc 
that is identical with the quantum Hamiltonian H ^ of Eq. (|29|) in the appropriate position- 
momentum representation, which may thus be diagonalised with the same rotation. The 

(2) 

equilibrium positions of Hsc are (+A X , —A y ) and (— A x , +A y ). 

The bounds on the existence of these classical effective Hamiltonians are exactly as in 
the quantum case - the excitation energies e_ and e_ of the decoupled modes remain 
real only on their respective sides of the critical coupling A c , which has the same value 
as in the quantum case. Clearly the semi-classical system is completely integrable in this 
thermodynamic limit. 



B. Equations of motion 

To analyse the behaviour of this semi-classical system for finite j, we form Hamilton's 
equations of motion from the derivatives of H sc [^J 

X = Px 

A [ZT xy 



y = Py 1 ~ 



2j V V 1 ~ V, 
p x = -u 2 x - 2Ay / tJcJo yy/l ~ V 



Py 



where as before, 



-uly - 2Xy/uJu^ xy/l-ri(l- — ^ ^ (68) 



It is not a priori obvious that this flow should preserve the condition set out in Eq. (|66|). 
However, we have demonstrated numerically that, providing we choose initial conditions 
that satisfy Eq. ([66)1 . then this condition is always satisfied. Although we have not shown 
this analytically, it can at least be seen to be plausible. Calculating rj = drj/dt = {H,rj} 
where {. . .} denote Poisson brackets, we find that 

A AuT 



V = — r\ — xpyy/1 - r], (70) 

J V Wo 

so that as r] approaches unity its rate of change approaches zero, implying that it is bound 
appropriately. 
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We now determine the fixed points of this flow at finite j by setting x = y = 0, p x = 
p y = 0. The simplest fixed point is given by x = y =p x = p y = 0, the co-ordinate origin. By 
calculating the Hessian stability matrix from the second derivatives of H, we see that this 
fixed point is only stable when 

\ + wg - V / (^_^) + i6A^ (l + l/(4j))} > 0. (71) 

i.e. when 

A < - Ac (72) 
y/l + WW) 

There are two other fixed points, both of which have p x = p y = 0, and with x and y given 
by 



2A 

Xq = ± 

UJ 



\ 





»o = T, Ht 1 + 77-Tf ■ (73) 



1 A 

These two quantities only remain real provided that 1 + — ^§ > 0, which corresponds to 



the condition 



A > , ^ , • (74) 



v /1 + V(4i) 

Provided that the above condition is fulfilled, the fixed points given by (+Xq, —yo) and 



-Xo,+yo) exist and are stable. So, below the coupling A c /yl + 1/ (4j), only one fixed 



point exists, which lies at the co-ordinate origin and is stable. Above A = A c / a/1 + 1/ (4j), 
this fixed point becomes unstable and two new stable fixed points appear at the coordinates 
(+xo, — yo) and (— xo,+yo). Note that these expressions give us the first correction to the 
location of the critical coupling in terms of a perturbation series in j. 

We can consider this semi-classical system as a particle moving in the two-dimensional, 
momentum-dependent potential 



U (x, y, p y ) = i (uj 2 x 2 + uly 2 ) + 2A y/uJuJ^ xy^jl ° V -. (75) 

Maps of this potential for different values of increasing coupling and for two different values 
of p y are shown in Fig. ED Firstly, note how greatly the value of p y affects the shape of the 
potential felt by the "particle". For example, above A c at A = 0.8, with p y = the potential 
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FIG. 13: The momentum-dependent potential U (x,y,p y ) at two different values of momentum 
p y = (left) and p y = 3 (right), for a series of couplings the same as in Fig. (|12j) . Note the 
difference in scales between the two plots. The Hamiltonian is resonant: u = ojq = 1, A c = 0.5 

bifurcates into two separate wells, whereas for p y = 3 it does not. Also note the similarity 
between the plot of U (x, y,p y ) for p y = and the plot of the wavefunction in Fig. H21 It is 
clear that the p y = potential largely determines the structure of the wavefunction at finite 
j, presumably because the location of the fixed points are determined with p y = 0. 

C. Chaos in the Semi-Classical model 

We numerically integrate Hamilton's equations of motion for the semi-classical system for 
a variety of different parameters and initial conditions. In order to analyse the trajectories 
resulting from these integrations, we use Poincare sections through the four- dimensional 
phase-space. As this system is Hamiltonian, the energy, 

E = -ju + - (u 2 x 2 + p 2 x -u + uly 2 +p 2 y - LO ) + 2X^/uJuJ^ xy-s/X - rj (76) 

is conserved, and thus we define our surface of section by p x = with p y being fixed by the 
energy E. We only record traversals for p y > 0. Poincare sections for illustrative parameter 
values are shown in Fig. 

At low A (A < 0.4 in Fig. I14p . the Poincare sections consist of a series of regular, periodic 
orbits. Approaching the critical coupling (A = 0.44,0.5 in Fig. I14jh we see a change in the 
character of the periodic orbits and also the emergence of a number of chaotic trajectories. 
Increasing the coupling further results in the break up of the remaining periodic orbits and 
the whole phase space becomes chaotic for couplings a little over the critical value (A = 0.6 

38 




FIG. 14: Poincare sections for the classical Dicke Model for a sequence of increasing couplings, 
with j = 5 and E = — 3. The Hamiltonian is resonant u = uo = 1;A C = 0.5 

in Fig. IT4"|) . This transition to chaos in the classical system mirrors very closely that seen 
in the quantum system, especially in the way that most of the change in the nature of the 
behaviour is centred about the critical coupling determined by the phase-transition. 

An interesting feature of this classical Hamiltonian is that the (re-)quantisation of this 
Hamiltonian is not unique. This is because the potential U(x,y,p y ) depends on the mo- 
mentum p y , a situation which may be compared to the quantisation of a Lagrangian for an 
electron in a magnetic field, where an extra 'rule' is required to obtain the correct quantisa- 
tion. We may requantise H sc by simply reversing the steps in Eqs f)61H65j) . However this is 
not the most obvious path as it involves the addition of extra imaginary p^-dependent terms 
which have canceled in the final Hamiltonian. Alternatively, one may simply requantise Eq. 
(Iri5j) as it stands, which results in the Hamiltonian 

/ Ub 

H' = ooo {b ] b - j) + uata + A (a f + a) {tf + b) J 1 - — , (77) 

which is clearly different to the original bosonic Hamiltonian of Eq. (|61|). This ambiguity 
disappears in the thermodynamic limit as here U (x, y,p y ) becomes momentum independent 
in this limit in both of the systems phases. 
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We note that the classical Hamiltonian 

H" = -ju + - (wV +PI-U + uly 2 + P 2 y - LU Q ) 

+ 2\^JU- X y^l-^Z^ (78) 

which is the same as the original Hamiltonian of Eq. (jHHj) but with p 2 removed from the 
square root, displays similar behaviour to that of the full Hamiltonian. The gain in simplicity 
in using this model, suggests that it would be an ideal test model for further exploration of 
the dynamics of this type of Hamiltonian constrained by a square-root. The behaviour of 
the Hamiltonian H" and fact that the p y = potential largely dominates dynamics of H c i 
and the structure of the wavefunction of the original DH, suggest that the requantisation 
route is not critical provided that j is not small. 
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VI. THE RWA AND INTEGRABILITY 



The DH in the RWA is given by 

Hrwa = u J z + uja)a + — ^= {a) J_ + aJ + ) . (79) 

It is in this form that the DH is generally studied and in which the thermodynamics of 
the phase transition were originally discussed [3, 2f|. In the RWA, the OPT occurs at 
a coupling twice that of the non- RWA critical value A* WA = 2A C = ^jUT Q Q This 
is simply a consequence of the fact that in the non-RWA DH there are four terms in the 
interaction, whereas here we only have two. As each term contributes to the mean-field, the 
critical coupling of the RWA is twice as big as the non-RWA one. 

In the RWA, the excitation number N of Eq. (jHj) becomes exactly conserved. This 
splits the total Hilbert space into an infinite number of sub-spaces, labeled by the excitation 
number n = 0, 1, 2, . . ., which in turn leads to level crossings and to a Poisson distribution 
for the P(S). The crossover between the RWA and non-RWA P (S) distributions has been 
studied by treating the non-RWA terms as a perturbation [34]], and it was found that as the 
strength of this perturbation is increased from zero to one, a standard crossover between 
Poissonian and Wigner-Dyson statistics is observed. 

Here we wish to report two observations concerning the difference between the RWA 
and non-RWA models. Firstly; a calculational issue that arises when considering the RWA 
system in the thermodynamic limit. We may derive effective Hamiltonians in each phase, 
by using the Holstein-Primakoff representation as before. In the normal phase, we obtain 

Hj£jy A = u Q b ] b + ua ] a + A (a ] b + b ] a) - joj - (80) 

The Bogoliubov transformations required to diagonalise this Hamiltonian are much simpler 
in terms of annihilation and creation operators than those for the non-RWA case. Specifi- 
cally, the RWA diagonalising transformations are 

a — > — c\ sin/5 + C2 cos/5; b — > c\ cos (3 + ci sin /3, (81) 

plus the Hermitian conjugate relations, where the rotation angle j3 is given by 

2A 

tan {2j3) = . (82) 

UJ — UJq 
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FIG. 15: The full energy schema of the (a) non-RWA and (b) RWA Dicke Hamiltonian for j=5. 
The Hamiltonian is resonant; uj = ujq = 1. 

The transformation for annihilation operators only involves annihilation operators, and 
the same with the creation operators. This is to be contrasted with the non-RWA transfor- 
mations, which transform any given operator into a linear combination of all four operators. 
Therefore, in the RWA it is very simple to find the diagonalising transformation in the second 
quantised representation, where as in the non-RWA case, this diagonalisation only becomes 
transparent when one considers the first quantised position-momentum representation of 
the operators. The converse of this statement is true; it is hard to find the diagonalising 
transformation in the RWA if one works in the position-momentum representation. We 
conjecture that this is a more general point than just applying here, and hope that this 
observation may be useful in other problems. 

Our second observation concerns the comparison of the energy spectra at finite j of the 
RWA and non-RWA Hamiltonians. Figure E3 shows two typical spectra, with coupling 
axes chosen for easy comparison. In terms of the appropriate critical coupling, the ground- 
state energy of the non-RWA spectrum is remarkably well approximated by the caustic of 
all the energy levels in the RWA spectrum that have negative slopes. As j increases this 
approximation becomes better, as the length of the individual line segments become shorter 
until, in the thermodynamic limit, the correspondence of the ground-states becomes exact 
and both excitation spectra become quasi-continuous. 
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VII. DISCUSSION 



We have presented a coherent and comprehensive picture of how the existence of a QPT 
in the thermodynamic limit plays a crucial role in determining chaotic properties in a model 
interacting system. The DH exhibits a change-over from quasi-integrability to chaos, and 
this transition is located by the precursors of the QPT around the critical point, A c . This 
statement applies equally well to the original quantum system and to the semi-classical 
counterpart derived from it. 

Our analysis of the DH in the thermodynamic limit consists of deriving an effective 
Hamiltonian to describe the system in each of its normal and super-radiant phases. For 
arbitrary coupling the system is described in terms of two decoupled modes, each of which is 
a collective photon-atom excitation, and it is the vanishing of the excitation energy associated 
with the photon-like mode that delimits the two phases. Our approach is particularly useful 
because we can calculate exactly any property of the system in the thermodynamic limit, 
by simply utilising the appropriate Bogoliubov transformations. 

This analysis reveals that the QPT breaks the symmetry associated with the parity op- 
erator II. In the normal phase, where the system in effectively unexcited, the wavefunctions 
of the system are invariant with respect to n. In the super-radiant phase however, this 
global symmetry is broken and two new local symmetries appear, each of which describes an 
isolated wavefunction lobe, and the spectrum is doubly degenerate. This symmetry breaking 
strictly only occurs in the thermodynamic limit, and at any finite j, these lobes are joined 
together in a total wavefunction that is II- invariant. That these two lobes are separated by 
a macroscopic amount, proportional to the square-root of the system size, means the onset 
of chaos is accompanied by the delocalisation of the wavefunctions and the appearance of 
macroscopic coherence in the system. 

Similar features occur in the three-dimensional Anderson model. This model of a disor- 
dered electron system exhibits a metal-insulator transition, in which the wavefunctions are 
localised for strong disorder and delocalised when the disorder is weak. Analysis 
of the level-statistics shows that the P (S) changes from Poissonian to Wigner-Dyson at 
the phase-transition point, which is determined by the magnitude of the random potential 
fluctuations. It is remarkable that our comparatively simple model should bear so many 
important features in common with complex disorder models, such as the Anderson model, 
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although one feature of such models that we have found no evidence of in the DH is the 
existence of a third universal P (S) distribution precisely at the critical coupling [3]. 

There are two different classical limits involved with the Dicke model, and by extension, 
models of similar nature involving quantum spins and boson fields. Firstly there is the limit 
of j — > oo, in which the length of the spin becomes macroscopic. The second is the limit 
h — > 0, which we have performed here when setting bosonic commutators equal to zero. 

These limits may be applied independently and in either order. If we apply the j — > oo 
limit first to the DH, we obtain the effective Hamiltonians H^ 1 ' 2 ' . Taking then h — > 
by setting the commutators of the collective modes to zero, we simply obtain iJ sc ' , the 
two classical effective Hamiltonians. Note that the integrability of H^ 1,2 ' makes this "de- 
quantisation" direct and unambiguous. Applying this limit in the other order means that 
starting with the DH in the Holstein-Primakoff representation we set h — > by setting the 
original field and atom bosonic commutators to zero, which results in our semi-classical 
Hamiltonian H sc . Subsequent taking of the j 1 — > oo limit results in if sc ' as above, showing 
that we obtain the same result independently of the order in which the limits are taken. 

fl 2) 

After both limits, the system described by if sc ' is "the classical" analogue of the DH, 
describing a macroscopic collection of atoms in terms of classical variables. This system is 
completely integrable, and there is no sign of chaos either in it or its quantised counterpart 

#(1,2). 

These results support the recent argument put forward by Ballentine concerning the 
existence of so-called "semi-quantum chaos" jfoj. Semi-quantum chaos is that which arises 
from the coupling of a quantum and a classical system, neither of which are by themselves 
chaotic. Ballentine studied a model of a massive particle of mass m interacting with a 
spin-half. By considering the semi-classical limit of m —>■ oo, the semi-quantum system of 
a quantum spin interacting with a classical particle was realised. He demonstrated that 
as m — > oo, the chaos in the system rapidly disappeared. Our results here may be seen 
as the complement to this system, where the mass is kept constant but the length of the 
pseudo-spin is taken to the classical limit j — ■> oo. Given the integrability of the DH is this 
limit, there is certainly no semi-quantum chaos in our system, which lends additional weight 
to Ballentine's claim that semi-quantum chaos does not exist. 

The question then arises what is the status of the two systems obtained by only taking 
one of the two limits. In the case of only taking the j — > oo limit, the answer is simple; 
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#(1,2) 

is a direct quantisation of Hsc and describes quantum fluctuations around classical 
mean-fields. More interesting is the status of H sc . We have shown here that its behaviour 
matches very closely that of the quantum DH, and that it has been derived in an almost 
canonical way, so its mathematical status as the semi-classical counterpart of the DH seems 
reasonably secure, but what the relevance of this model to the physical system is less obvious. 

The nature of the H — > limit suggests that this model might be useful in describing the 
model when there are a few atoms (10-20) present, and almost-classical fields, i.e. coherent 
states, are applied. Under these circumstances the original DH and semi-classical model H sc 
might be fruitfully compared. 
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APPENDIX A: BOGOLIUBOV TRANSFORMATION 



1. Normal phase 



The two sets of bosons, {a, b} and {ci, c 2 }, may be expressed in terms of one another as 
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The angle 7^ is rotation angle of the coordinate system which eliminates the interaction in 
the position representation, and is given by 



tan (2 7 W) = 



(A3) 
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2. Super-radiant phase 



,t 



The analogous Bogoliubov transfomations is the super-radiant phase are 
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where the angle 7^ is given by 



tan (2 7 < 2 >) = 



and where 
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APPENDIX B: SQUEEZING VARIANCES 

The preceeding Bogoliubov transformations may be used to derive exact expressions for 
the squeezing variances of the ground-state wavefunction in the thermodynamic limit as 
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discussed in section IIHI In the normal phase they are given by 
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whereas in the super-radiant phase we find 
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These results are plotted in the main body of the text. 
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